In this lab, we’ll explore the basics of working with raster data, including attribute, spatial, and geometry operations. This lab follows chapters 3, 4, and 5 of Geocomputation with R by Robin Lovelace.

Libraries Included

library(terra)
## Warning: package 'terra' was built under R version 4.1.2
## terra 1.6.17
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.1.2
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:terra':
## 
##     intersect, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(spData)
## Warning: package 'spData' was built under R version 4.1.2
library(spDataLarge)
library(tmap)
## Warning: package 'tmap' was built under R version 4.1.2
library(geodata)
## Warning: package 'geodata' was built under R version 4.1.2

Raster Data

Geometry Includes the following components

Attributes contain the following components

Subsetting for Raster Data

Packages in R’s Spatial Ecosystem: raster (core package for rasters), sf (new and more popular), sp (original package), stars (can handle both vector and rasters), terra (raster manipulation)

Manipulating Raster Objects

Loading in a raster object

Raster data represents continuous surfaces, as opposed to the discrete features represented in the vector data model. Here we’ll learn how to create raster data objects from scratch and how to do basic data manipulations.

Let’s create a SpatRaster object using a digitial elevation model for Zion National Park.

raster_filepath <- system.file("raster/srtm.tif", package = "spDataLarge")
my_rast <- rast(raster_filepath)

class(my_rast) # test class of raster object
## [1] "SpatRaster"
## attr(,"package")
## [1] "terra"
#checking what the raster is and giving summary information
my_rast
## class       : SpatRaster 
## dimensions  : 457, 465, 1  (nrow, ncol, nlyr)
## resolution  : 0.0008333333, 0.0008333333  (x, y)
## extent      : -113.2396, -112.8521, 37.13208, 37.51292  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : srtm.tif 
## name        : srtm 
## min value   : 1024 
## max value   : 2892
#As you can see this is a new class: SpatRaster, with specific dimensions included
#Dimensions show that there are 457 rows, 465 columns, and nlyr indicates how many layers are included (in this case nlyr = 1)
#Resolution indicates the units that the raster resolution is in

Let’s quickly polot the raster

#Using plot
plot(my_rast)

#Using tmap

tm_shape(my_rast) +
  tm_raster()

Creating a raster data from scratch

We can also create rasters from scratch using the rast() function. Here we create 36 cells centerd around (0,0). By default the CRS is set to WGS84, but we could change this with the “crs” argument. Because we are working in WGS84, the resolution is in units of degrees. rast() fills the values of the cells row-wise starting in the upper left corner.

elev <- rast(nrows = 6, ncols = 6, resolution = 0.5, #setting the columns and rows and resolution
             xmin = -1.5, xmax = 1.5, ymin = -1.5, ymax = 1.5, #setting the extent/bounding box
             vals = 1:36) #creating a matrix of 6x6 to create 36 grid cells filled 1-36 with 1 step for each pixel
#You can define the CRS as well but we are just keeping it at WGS84 for now

#Checking our output raster
elev
## class       : SpatRaster 
## dimensions  : 6, 6, 1  (nrow, ncol, nlyr)
## resolution  : 0.5, 0.5  (x, y)
## extent      : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source      : memory 
## name        : lyr.1 
## min value   :     1 
## max value   :    36
#Plot output raster
plot(elev)

Pulling in Landsat Data

The SpatRaster class can also handle multiple layers.

#pulling in landsat
multi_raster_file <- system.file("raster/landsat.tif", package = "spDataLarge")

#turning this into a raster
multi_rast <- rast(multi_raster_file)

#checking what this raster is
multi_rast
## class       : SpatRaster 
## dimensions  : 1428, 1128, 4  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 301905, 335745, 4111245, 4154085  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 12N (EPSG:32612) 
## source      : landsat.tif 
## names       : landsat_1, landsat_2, landsat_3, landsat_4 
## min values  :      7550,      6404,      5678,      5252 
## max values  :     19071,     22051,     25780,     31961
#names indicate the different landsat bands, there are 4 layers included

#Checking the number of layers included in the raster object 
nlyr(multi_rast)
## [1] 4

Subsetting a band from Landsat

  • We can subset layers using either the layer number or name

  • We can combine SpatRaster objects into one, using the c function

#pulling out band 3
multi_rast3 <- subset(multi_rast, 3)
nlyr(multi_rast3)
## [1] 1
#pulling out band 4
multi_rast4 <- subset(multi_rast, "landsat_4")
nlyr(multi_rast4)
## [1] 1
multi_rast34 <- c(multi_rast3, multi_rast4)

Plotting full multi-layer raster

plot(multi_rast)

Character Vectors and assigning locations with soil types

grain_order <- c("clay", "silt", "sand")
grain_char <- sample(grain_order, 36, replace = TRUE)
grain_fact <- factor(grain_char, levels = grain_order) #this is just 36 randomly selected silt, sand, clay (not a raster yet)

#Turning it into a raster 

grain <- rast(nrows = 6, ncols = 6, resolution = 0.5, #establishing the rows, columns and resolution
              xmin = -1.5, xmax = 1.5, ymin = -1.5, ymax = 1.5, #setting the extent/bounding box
              vals = grain_fact) #give it values from the grain factors we established before
#plot it
plot(grain)

Raster Subsetting

We can index rasters using, row-column indexing, cell IDs, coordinates, other spatial objects.

elev[1,1] #finding the value in row 1, column 1
##   lyr.1
## 1     1
elev[1] #finding the value in cell ID 1 (only putting in 1 value indicates cell id rather than rows and columns)
##   lyr.1
## 1     1
elev[36] #finding the value of cell 36
##   lyr.1
## 1    36

Creating multi-layer raster

two_layer <- c(grain, elev)
two_layer[1]
##   lyr.1 lyr.1.1
## 1  sand       1
#Updating one layer and replacing with new value 
elev[1,1] = 0 #updating the raster to manually override and replace the value that was in column 1,1 with 0
elev[1]
##   lyr.1
## 1     0
#updating two layers and replacing with new value
# Replacing values in multi-layer rasters requires a matrix with as many columns as layers and rows as replaceable cells.
two_layer[1] <- cbind(c(1), c(4))
two_layer[1]
##   lyr.1 lyr.1.1
## 1  silt       4

Summary Statistics for Raster Objects

We can get info on raster values just by typing the name or using the summary function.

#We can get global summaries, such as standard deviation.
summary(elev)
##      lyr.1      
##  Min.   : 0.00  
##  1st Qu.: 9.75  
##  Median :18.50  
##  Mean   :18.47  
##  3rd Qu.:27.25  
##  Max.   :36.00
global(elev, sd) 
##             sd
## lyr.1 10.58432
#Or we can use freq() to get the counts with categories.
freq(grain) #this tells us the number of grid cells for each value (count of number of pixels for each value 
##   layer value count
## 1     1  clay    10
## 2     1  silt    10
## 3     1  sand    16
hist(elev) #looking at a histogram of the values within the raster dataset

Spatial Subsetting

We can move from subsetting based on specific cell IDs to extract info based on spatial objects.

To use coordinates for subsetting, we can “translate” coordinates into a cell ID with the terra function cellFromXY() or terra::extract().

id <- cellFromXY(elev, xy = matrix(c(0.1, 0.1), ncol = 2))
elev[id]
##   lyr.1
## 1    16
#This is the same as
terra::extract(elev, matrix(c(0.1,0.1), ncol = 2))
##   lyr.1
## 1    16

Raster objects can also subset with another raster object. Here we extract the values of our elevation raster that fall within the extent of a masking raster.

clip <- rast(xmin = 0.9, xmax = 1.8, ymin = -0.45, ymax = 0.45,
             resolution = 0.3, vals = rep(1,9))

elev[clip] #extracting values from elev that fall within the clip mask that we created above 
##      lyr.1
## [1,]    18
## [2,]    24
terra::extract(elev, ext(clip))
##   lyr.1
## 1    18
## 2    24

In the previous example, we just got the values back. In some cases, we might want the output to be the raster cells themselves. We can do this use the “[” operator and setting “drop = FALSE”

This example returns the first 2 cells of the first row of the “elev” raster.

elev[1:2, drop = FALSE]
## class       : SpatRaster 
## dimensions  : 1, 2, 1  (nrow, ncol, nlyr)
## resolution  : 0.5, 0.5  (x, y)
## extent      : -1.5, -0.5, 1, 1.5  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source      : memory 
## name        : lyr.1 
## min value   :     0 
## max value   :     2

Masking one raster with another

Another common use of spatial subsetting is when we use one raster with the same extent and resolution to mask the another. In this case, the masking raster needs to be composed of logicals or NAs.

# create raster mask of the same resolution and extent
rmask <- elev

#randomly replace values with NA and TRUE to use as a mask
values(rmask) <- sample(c(NA,TRUE), 36, replace = TRUE)

#spatial subsetting 
elev[rmask, drop = FALSE] # with the [] operator, you can mask one raster with the other, in this case masking elevation with rmask 
## class       : SpatRaster 
## dimensions  : 6, 6, 1  (nrow, ncol, nlyr)
## resolution  : 0.5, 0.5  (x, y)
## extent      : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source      : memory 
## name        : lyr.1 
## min value   :     0 
## max value   :    36
#this is where the actual masking is being done
mask(elev, rmask)
## class       : SpatRaster 
## dimensions  : 6, 6, 1  (nrow, ncol, nlyr)
## resolution  : 0.5, 0.5  (x, y)
## extent      : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source      : memory 
## name        : lyr.1 
## min value   :     2 
## max value   :    35

Map Algebra

Here we define map algebra as the set of operations that modify or summarize raster cell values with reference to surrounding cells, zones, or statistical functions that apply to every cell.

Local Operations

Local operations are computed on each cell individually. We can use ordinary arithemetic or logical statements

#adding two rasters together
elev + elev
## class       : SpatRaster 
## dimensions  : 6, 6, 1  (nrow, ncol, nlyr)
## resolution  : 0.5, 0.5  (x, y)
## extent      : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source      : memory 
## name        : lyr.1 
## min value   :     0 
## max value   :    72
#squaring a raster
elev^2
## class       : SpatRaster 
## dimensions  : 6, 6, 1  (nrow, ncol, nlyr)
## resolution  : 0.5, 0.5  (x, y)
## extent      : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source      : memory 
## name        : lyr.1 
## min value   :     0 
## max value   :  1296
#taking the log
log(elev)
## class       : SpatRaster 
## dimensions  : 6, 6, 1  (nrow, ncol, nlyr)
## resolution  : 0.5, 0.5  (x, y)
## extent      : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source      : memory 
## name        :    lyr.1 
## min value   :     - ?  
## max value   : 3.583519
elev > 5
## class       : SpatRaster 
## dimensions  : 6, 6, 1  (nrow, ncol, nlyr)
## resolution  : 0.5, 0.5  (x, y)
## extent      : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source      : memory 
## name        : lyr.1 
## min value   : FALSE 
## max value   :  TRUE

Classifying intervals of values into groups

We can also classify intervals of values into groups. For example, we could classify a DEM into low, middle, and high elevation cells

#This is creating a reclassification matrix

rcl <- matrix(c(0,12,1,12,24,2,24,36,3), ncol = 3, byrow = TRUE)
#this is saying that there are three columns so the first column is going to have a 0, then 12, then 1, and then it will go to the next row

rcl
##      [,1] [,2] [,3]
## [1,]    0   12    1
## [2,]   12   24    2
## [3,]   24   36    3
#using this maxrix to then reclassify our elevation matrix
recl <- classify(elev, rcl = rcl) #we named rcl as our reclassification matrix, but if the reclassification matrix was named something else then it would be rcl = "name"
recl
## class       : SpatRaster 
## dimensions  : 6, 6, 1  (nrow, ncol, nlyr)
## resolution  : 0.5, 0.5  (x, y)
## extent      : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source      : memory 
## name        : lyr.1 
## min value   :     0 
## max value   :     3

For more efficient processing, we can use a set of map algebra functions.

We can use the lapp() function to compute the Normalized Difference Vegetation Index (NDVI). Let’s calculate NDVI for Zion National Park using multispectral satellite data.

multi_raster_file <- system.file("raster/landsat.tif", package = "spDataLarge")
multi_rast <- rast(multi_raster_file)

We need to define a function to calculate NDVI

#creating an ndvi function with individual layers within the image stack as the inputs
ndvi_fun = function(nir,red){
  (nir-red) / (nir + red)
}

So now we can use lapp() to calculate NDVI in each raster cell. To do so, we just need the NIR and red bands. As a reminder, lapp() applies a function to each cell using layers as arguments

#this is using lapp to apply the ndvi function to each cell
ndvi_rast <- lapp(multi_rast[[c(4,3)]], fun = ndvi_fun)

tm_shape(ndvi_rast)+
  tm_raster()
## stars object downsampled to 888 by 1125 cells. See tm_shape manual (argument raster.downsample)
## Variable(s) "NA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

Focal Operations

Local operations operate on one cell, though from multiple layers. Focal operations take into account a central (focal) cell and its neighbors. The neighborhood (or kernel, moving window, filter) can take any size or shape. A focal operation applies an aggregation function to all cells in the neighborhood and updates the value of the central cell before moving on to the next central cell

We can use the focal() function to perform spatial filtering. We define the size, shape, and weights of the moving window using a matrix. Here we find the minimum.

elev <- rast(system.file("raster/elev.tif", package = "spData"))

#The function being used here is the minimum function
r_focal <- focal(elev, w = matrix(1, nrow = 3, ncol = 3), fun = min)

plot(elev)

plot(r_focal)

Zonal Operations

Similar to focal operations, zonal operations apply an aggregation function to multiple cells. However, instead of applying operations to neighbors, zonal operations aggregate based on “zones”. Zones can are defined using a categorical raster and do not necessarily have to be neighbors

For example, we could find the average elevation for different soil grain sizes.

#zones where the statistic "mean" is being applied to elevation by first grouping all the soil types together (group the soil types together first and then take the mean elevation for each of those groups so we have average elevation for each type of soil)
zonal(elev, grain, fun = "mean")
##   lyr.1    elev
## 1  clay 22.5000
## 2  silt 21.4000
## 3  sand 14.1875

Merging Rasters

aut = geodata::elevation_30s(country = "AUT", path = tempdir())
ch = geodata::elevation_30s(country = "CHE", path = tempdir())
aut_ch = merge(aut, ch) #merging these two raster datasets together

Geometric Operations

elev = rast(system.file("raster/elev.tif", package = "spData"))
elev_2 = extend(elev, c(1, 2)) # add one row and two columns

plot(elev)

plot(elev_2)

Performing algebraic operations on objects with different extents doesn’t work.

#this doesn't work because the two extents do not match
#elev + elev_2

We can align the extent of the 2 rasters using the extend() function. Here we extend the elev object to the extent of elev_2 by adding NAs.

elev_4 <- extend(elev, elev_2)

the origin function returns the coordinates of the cell corner closes to the coordinates (0,0). We can also manually change the origin.

#manually changing the orgin 
origin(elev_4)
## [1] 0 0
origin(elev_4) <- c(0.25, 0.25)
origin(elev_4)
## [1] 0.25 0.25

Aggregation and Disaggregation

Faster datasets can also differ in their resolution to match resolutions we can decrease the resolution by aggregating or increase the resolution by disaggregating.

Let’s start by changing the resolution of a DEM by a factor of 5, by taking the mean.

#Aggregating

dem <- rast(system.file("raster/dem.tif", package = "spDataLarge"))
dem_agg <-  aggregate(dem, fact = 5, fun = mean)

plot(dem)

plot(dem_agg)

Increasing the Resolution

We have some choices when increasing the resolution. Here we try the bilinear method.

#increasing the resolution using method = "bilinear"
dem_disagg <- disagg(dem_agg, fact = 5, method = "bilinear")
identical(dem, dem_disagg)
## [1] FALSE
plot(dem_disagg)

Resampling

Aggregation/disaggregation work when both rasters have the same origins. What do we do in the case where we have two or more rasters with different origins and resolutions? Resampling computes values for new pixel locations based on custom resolutions and origins.

In most cases, the target raster would be an object you are already working with, but here we define a target raster.

target_rast <- rast(xmin = 794600, xmax = 798200,
                   ymin = 8931800, ymax = 8935400,
                   resolution = 150, crs = "EPSG:32717")

dem_resampl <- resample(dem, y = target_rast, method = "bilinear")

plot(dem)

plot(dem_resampl)